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A method for calculating the pressure tensor in const ant- volume Monte Carlo simulations of con- 
vex bodies is presented. In contrast to other approaches, the method requires only an isotropic 
scaling of the simulation box, and the counting of simple geometric quantities characterizing over- 
lapping pairs. Non-sphericity presents no special difficulties. The result is expressed as a sum of 
pairwise contributions, and can therefore be used to compute pressure tensor profiles in a conven- 
tional way. 
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I. INTRODUCTION 

The purpose of this brief paper is to explain in 
detail a procedure to calculate the pressure tensor 
in constant-volume computer simulations of hard 
convex bodies. The term "hard" is used to de- 
note particles whose interaction energy is infinite 
when they overlap, and zero when they do not; 
this paper will also consider "soft" convex bodies, 
meaning particles whose interaction energy takes 
a finite, positive, constant value e when they over- 
lap, and zero when they do not. The restriction 
to convex bodies is a simplification to ensure that 
an overlap between a pair of molecules may result 
when the distance between the molecular centres is 
reduced, but not when it is increased, while keep- 
ing all orientations fixed; extending the results to 
non-convex bodies is relatively straightforward but 
requires care. Calculation of the components of 
the pressure tensor, as opposed to the scalar pres- 
sure, is an important route to the surface tension 
of interfaces. 

The motivation for this work is the recent pa- 
per by Gloor et al.- which, as well as providing 
a near-comprehensive review of simulation tech- 
niques in this area, proposes a method for surface 
tension calculation based on perturbations in the 
cross-sectional area of a system containing the pla- 
nar interface of interest. It is claimed that a key 
advantage of this "test-area" method for discon- 
tinuous potentials is the avoidance of the determi- 
nation of delta-functions "which are impractical to 
evaluate, particularly in the case of non-spherical 
molecules" . The present paper demonstrates that 
this is not the case: that the pressure tensor com- 



ponents for such systems, and hence the surface 
tension, may be calculated by a simple extension 
of the standard approach of counting overlaps pro- 
duced by an isotropic volume scaling. In fact, this 
method has already been used to determine the 
surface tension of the isotropic-nematic interface 
of hard ellipsoids, 2 although the full details of the 
method were not explained in that paper. 

Microscopic expressions for the pressure in a 
fluid are to be found in the standard referencesi&i 
Separate P = P id + P cx where P id = pk B T is 
the ideal gas contribution, and focus on the ex- 
cess part, P cx . The usual mechanical route to the 
pressure starts with the volume derivative of the 
excess Helmholtz free energy F cx 
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where T is the temperature and k B Boltzmann's 
constant. The excess partition function in the 
canonical ensemble, Q cx , is written 
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where U is the total potential energy. The iV parti- 
cle centre-of-mass coordinates are denoted rj , and 
scaled coordinates s j are defined by rj — Lsj , as- 
suming a cubic box of side L and volume V — L 3 . 
Particles may have orientational degrees of free- 
dom; for example the orientation of a uniaxial 
molecule is specified by a unit vector Uj directed 
along the symmetry axis. However, for simplicity 
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here, and in most of what follows, the orientational 
dependence and the integrations over orientational 
degrees of freedom are not explicitly written. In 
the following, also for simplicity, pairwise-additive 



potentials U = ^ 



i<j Ui J 



with Uij = u(rij,Ui,iij), 



where r y = r; — r ? - , are assumed. For continuous 
potentials 



the volume differentiation within 
the integral of eqn (J2J) is readily carried out yield- 
ing the well-known expression 
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Any orientations (for example Ui, Uj, and = 
r ij/ r ij) are held fixed in taking the derivative in 
the right-most expression in eqn J^J. 
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All pairs are now equivalent in this expression. For 
hard particles, each exponential is either (sig- 
nalling a newly-generated overlap) or 1 (when the 
volume scaling does not cause an overlap), and 
the average has the interpretation of a probability 
that any given pair will not overlap. For soft par- 
ticles, every newly-generated overlap contributes 
e -e/k B T ^ w jjjj e p a j rs whose overlap status remains 
unchanged contribute 1. In either case, it is con- 
venient to define 



II. PRESSURE FOR DISCONTINUOUS 
POTENTIALS 

For discontinuous potentials, dQ ex /dV may be 
estimated by finite differences. To highlight the 
connection with the method for determining the 
full pressure tensor, a derivation is presented here; 
this follows Eppenga and Frenkel-, Perram et al^, 
and the basic ideas go back to Vieillard-Baroni, 
although the detailed approach and notation used 
here are different. Consider a volume change V — > 
V - AV. 



P cx _ 1 dQ cx 
fc^T ~ Q cx dV 
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The ensemble average is conducted in the unper- 
turbed system and AU = Ei<j ^ u ij IS a sum °f 
energy changes for each pair, arising from overlaps 
which are freshly generated by the volume change. 
Hence the Boltzmann factor is expressed as a prod- 
uct of terms. The expression is simplified by the 
assumption that the different pair overlaps are un- 
corrected. This is reasonable at small AV, since 
the number of overlaps may be made vanishingly 
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■poveriap will be loosely termed an 



overlap prob- 
ability" ; for small AV, it is small compared with 



1, simply because most pairs will have Aui 



0. 



The indices i and j on -p° vcrla P are actually redun- 
dant. The product of terms is expanded to first 
order in the overlap probabilities: 
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For hard particles, $° voria P = 7V ovorla P, the num- 
ber of overlaps in a configuration generated by the 
volume scaling, arising from summing the proba- 
bilities over all pairs. In the more general case of 
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soft particles 

^overlap = V y = V 1 - Q-^ij/knT 

>(l-e 



re-emphasizing that 7V ovcrlap counts the number of 
new overlaps, since the original configuration may 
already contain overlaps if e is finite, and these do 
not contribute. 

Suppose now that the volume reduction comes 
from scaling the box lengths by a factor (1 — e), i.e. 
L -> L- AL where AL/L = e, then AV/V « 3e 
and 
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Thus the pressure may be estimated by essentially 
counting the overlaps that result from a small 
isotropic volume scaling. The same expression (to 
leading order in e) results if the overlaps result 
from scaling all linear particle dimensions by the 
factor 1 + e instead of reducing the volume. 

The expression of Perram et al^ for convex 
spheroids is equivalent. It may be written 
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with overlap function F^ defined such that 

Fij > 1, for non-overlapping i & j, 
F^ = 1, for i & j in contact, 
F^ < 1, for overlapping i & j. 

The sum in eqn Q counts all distinct pairs for 
which 1 < < (1 + e) 2 w 1 + 2e. The value F tJ 
within the sum can be replaced by 1. The overlap 
function defined by Perram et al.— has the scal- 
ing behaviour Fij — r 2 fij(fij,Ui,Uj) so counting 
pairs for which 1 < Fij <~ (1 + e) 2 is equivalent 
to counting the overlaps generated by scaling all 
particle linear dimensions by 1 + e as above. 




FIG. 1: Schematic geometry of two convex particles 
in contact (left) and near contact (right). The vectors 
£i and as well as the particle orientations, are kept 
fixed as the centre-centre vector rij is scaled uniformly. 
The separation distance is defined between the two 
tangent planes. After Ref. 0- 



III. PRESSURE TENSOR FOR 
DISCONTINUOUS POTENTIALS 



Turning now to the pressure tensor, the ideal 
be expressed^ 



part is P^a/ksT = p5 a p and the excess part may 
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where ig) represents an outer product. To estimate 
this tensor in the case of discontinuous potentials, 
a more detailed consideration of pair geometry is 
required. When the particles are in contact, there 
is a common tangent plane at the point of contact; 
define a unit vector Sij normal to this plane, from 
j to i, as shown in Fig.^ For small displacements 
of the centre-centre vector ry, keeping all angles, 
and the vectors £j and £j of the contact points rel- 
ative to the molecular centres, fixed, the contact 
points and their associated tangent planes sepa- 
rate. Define the separation vector £y , and the per- 
pendicular separation between the tangent planes 

= iij ■ £ij , as shown in the figure. Then 
takes positive values if the particles do not overlap, 
and negative values if they do overlap. Hence the 
step function 9(sij) takes the value (overlap) or 
1 (no overlap), and a standard trick^-' 11 ! 12 may 
be used to estimate the gradient of the potential. 
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Write 



-Uij/k B T _ 



( Sij ) + (1 - 9( Sij ))e 



-e/k B T 



(the second term is zero for hard particles). Take 
the gradient of both sides and rearrange 
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Note that the delta function may be taken to act 
at a separation infinitcsimally outside contact, and 
so the exponential e +Uij / kBT may be replaced by 

unity. Also, [ — — ] = §ij, so the expression sim- 



plifies: 



dri 



This is a straightforward generalization of the ex- 
pression of Boublik 9 to include the case of soft 
particles. The connection to the pair distribution 
functions of hard convex body fluids at contact is 
set out elsewhere^ For the present purposes, the 
delta function is approximated^ by the term in 
braces below: 

1 du[Tjj) _ 

k B T dra ~ SlJ 



x < lim ■ 
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Arij with Arij 



As, 



If the change results from a scaling r, 
trij, then As = Si 
eSij ■ nj . The remaining terms may be 
recognized as (f>ij, defined in eqn Q. The result 
becomes 
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This is the main result: pressure tensor compo- 
nents in constant-volume Monte Carlo simulations 
of general convex bodies can be calculated al- 
most as simply as the scalar pressure, by sum- 
ming over pairs which would overlap as a result 
of an isotropic scaling of coordinates. The term 



4>ij acts to select such pairs (and for soft parti- 
cles, it also incorporates the appropriate energy- 
dependent scale factor); the only additional re- 
quirement is to compute the surface normal for 
the near-contact pairs, which is almost always 
straightforwardly done. The above derivation fol- 
lows closely that of Trokhymchuk and Alejandro^, 
although that paper explicitly considered only the 
spherically symmetric case. Taking one third of 
the trace of eqn JHJ regenerates eqn © . 

This expression is consistent with the form of 
pressure tensor calculated in collisional molecular 
dynamics simulations 



colls 



which is a sum over collisions occurring in time t, 
where the colliding particles are labelled i and j, 
and the impulsive momentum transfer of magni- 
tude Apij is directed along the common surface 
normal sy. For continuous potentials, the analo- 
gous formula is proportional to ® where 
is the pair force. Indeed, eqn JHJ may be written 
down almost immediately, knowing that pair con- 
tributions must be proportional to 
that the isotropic pressure is given by eqn (jSJ 



TV,-, and 



To evaluate the surface tension 7 of a planar in- 
terface normal to the z direction, the difference in 
components AP(z) = P zz - \{P xx {z) + P yy {z)) 
is calculated as a function of z and integrated 
across the interface. The integrals may be com- 
puted implicitly, yielding 7 as a single number, 
but the above formalism makes explicit the con- 
tributions of each pair to the total for each com- 
ponent of P, and therefore the definition of a local 
pressure tensor P(z) via the Irving-Kirkwood or 
Harasima conventioni^i£ii£ is handled in exactly 
the same way as for continuous potentials. The in- 
tegration of AP{z) is then performed numerically. 
As emphasized by Holcomb et alJ£ and Trokhym- 
chuk and Alejandre^, monitoring pressure tensor 
profiles in inhomogeneous systems is an important 
test of equilibrium, as well as providing useful in- 
formation about interface structure. This aspect 
is missing from the test-area method as described 
in Ref. although it could probably be incorpo- 
rated. 
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IV. CONCLUSIONS 

This paper has demonstrated that all compo- 
nents of the pressure tensor may be straight- 
forwardly computed in constant-volume simula- 
tions by counting the pair overlaps resulting from 
isotropic scaling of the simulation box, and includ- 
ing a simple geometrical tensor formed from the 
surface normal at contact and the centre-centre 
vector. Non-spherical particles with discontinuous 
interaction potentials present no great difficulties 
and special anisotropic scaling moves are not re- 
quired. The method makes it easy to resolve the 
pressure tensor profile, for use in surface tension 
calculations. 

No explicit comparisons of techniques are per- 
formed here, but it should be noted that equation 
(JSJ has been used to calculate pressure tensor pro- 
files, and the surface tension, of the equilibrium 
isotropic-nematic interface in the hard ellipsoid 
fluid, 2 and comparisons made there with results 
for a closely-related continuous potential model. 
There is a trade-off between minimizing statisti- 
cal errors (favoured by large numbers of overlaps) 
and systematic errors (favoured by small volume 



changes), in the choice of the scaling parameter e, 
so usually it will be necessary to examine several 
values. As has been emphasized, the method is es- 
sentially equivalent to a numerical estimate of the 
(orientationally-averaged) pair distribution func- 
tion at contact, and so the same caveats apply as 
to the calculation of that quantity^ near contact, 
the function may vary quite sharply, so the extrap- 
olation to small e may require care. These obser- 
vations apply equally to the test-area method, and 
all other finite perturbation methods of this kind. 
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